Streams
3.1 Import AOI
library(sf)
aoi_00 = geodata::gadm(
country="CAN", level=0,
path="./assets/SHP/") |>
sf::st_as_sf()
aoi_01 = geodata::gadm(
country="CAN", level=1,
path="../assets/SHP/") |>
sf::st_as_sf()
aoi_02 = geodata::gadm(
country="CAN", level=2,
path="../assets/SHP/") |>
sf::st_as_sf()
aoi_03 = geodata::gadm(
country="CAN", level=3,
path="../assets/SHP/") |>
sf::st_as_sf()
# ------------ #
# BC Albers #
# ------------ #
crs_master = "EPSG:3005"
sf::st_transform(aoi_00, crs_master)
sf::st_transform(aoi_01, crs_master)
sf::st_transform(aoi_02, crs_master)
sf::st_transform(aoi_03, crs_master)
sf::st_write(aoi_00, "../assets/SHP/aoi_00.shp", delete_dsn = T)
sf::st_write(aoi_01, "../assets/SHP/aoi_01.shp", delete_dsn = T)
sf::st_write(aoi_02, "../assets/SHP/aoi_02.shp", delete_dsn = T)
sf::st_write(aoi_03, "../assets/SHP/aoi_03.shp", delete_dsn = T)3.2 Extract Shoreline
library(osmdata)
library(rmapshaper)
bc_province = sf::st_read("../assets/SHP/aoi_01.shp", quiet=T) |>
sf::st_make_valid() |> dplyr::filter(NAME_1 == "British Columbia")
bc_districts = sf::st_read("../assets/SHP/aoi_02.shp", quiet=T)|>
sf::st_make_valid() |> dplyr::filter(NAME_1 == "British Columbia")
bc_counties = sf::st_read("../assets/SHP/aoi_03.shp", quiet=T)|>
sf::st_make_valid() |> dplyr::filter(NAME_1 == "British Columbia")
# --------------------------- #
# Note: BBOX needs CRS forcing #
island_extent = sf::st_bbox(c(
xmin = -129.5, ymin = 48.1,
xmax = -123.0, ymax = 51.0),
crs = st_crs(4326)) |>
sf::st_as_sfc()
island_extent_sf = sf::st_as_sf(island_extent)
shoreline_input <- osmdata::opq(bbox = sf::st_bbox(island_extent_sf)) |>
osmdata::add_osm_feature(key = "natural", value = "coastline") |>
osmdata::osmdata_sf() |> {\(x) x$osm_lines}() |>
sf::st_intersection(island_extent_sf) |>
sf::st_union() |>
sf::st_polygonize() |>
sf::st_collection_extract("POLYGON")
# `rmapshaper::ms_clip()` requires perimiter polygon shell only
shoreline_shell <- shoreline_input %>%
sf::st_union() %>%
sf::st_as_sf()|>
sf::st_transform(crs_master)
island_shoreline <- bc_districts %>%
sf::st_filter(shoreline_input, .predicate = st_intersects) %>%
rmapshaper::ms_dissolve() %>%
rmapshaper::ms_clip(shoreline_shell)|>
sf::st_transform(crs_master)
island_districts <- bc_districts %>%
sf::st_filter(shoreline_input, .predicate = st_intersects) %>%
rmapshaper::ms_simplify(keep = 0.95) %>%
rmapshaper::ms_clip(shoreline_shell)|>
sf::st_transform(crs_master)
island_counties <- bc_counties %>%
sf::st_filter(shoreline_input, .predicate = st_intersects) %>%
rmapshaper::ms_simplify(keep = 0.95) %>%
rmapshaper::ms_clip(shoreline_shell)|>
sf::st_transform(crs_master)
sf::st_write(island_districts, "../assets/SHP/island_districts.shp", delete_dsn=T)
sf::st_write(island_counties, "../assets/SHP/island_counties.shp", delete_dsn=T)
sf::st_write(island_shoreline, "../assets/SHP/island_shoreline.shp", delete_dsn=T)
sf::st_write(shoreline_input, "../assets/SHP/shoreline_input.shp", delete_dsn=T)
sf::st_write(shoreline_shell, "../assets/SHP/shoreline_shell.shp", delete_dsn=T)
tmap::tmap_mode("view")
tmap::tm_shape(island_counties) + tm_borders(lwd=0) +
tmap::tm_shape(island_shoreline) + tm_borders(col = "black", lwd=4) +
tmap::tm_shape(island_districts) + tm_borders(col = "purple", lwd=1.5) +
tmap::tm_shape(island_counties) + tm_borders(col = "yellow", lwd=0.5) +
tmap::tm_basemap("Esri.WorldImagery") -> tm01_live
tm01_live